Skewness Test (D’Agostino)#
The skewness test checks whether a sample shows statistically significant asymmetry.
It is commonly used as a diagnostic for:
symmetry assumptions (many parametric procedures implicitly assume roughly symmetric errors)
whether a transformation (log/sqrt/Box–Cox) might be useful
as one component of normality diagnostics (but it is not a full normality test by itself)
Learning goals#
By the end you should be able to:
explain what skewness measures (right/left tail intuition)
state the hypotheses behind the skewness test
implement D’Agostino’s skewness test with NumPy only
interpret the test statistic and p-value (including one-sided alternatives)
understand common pitfalls (outliers, sample size, “significant but tiny” effects)
import numpy as np
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=6, suppress=True)
rng = np.random.default_rng(42)
Prerequisites#
Sample mean and central moments
Z-scores and p-values
The standard normal distribution
If any of these are fuzzy, you can still follow the notebook: the code is written to make each step explicit.
1) Intuition: what skewness measures#
Skewness is a signed measure of asymmetry:
Positive skewness → a longer right tail (“right-skewed”)
Negative skewness → a longer left tail (“left-skewed”)
A common (moment-based) sample skewness is
[ ext{skew}(x) = rac{m_3}{m_2^{3/2}}, \qquad m_k = rac{1}{n}\sum_{i=1}^n (x_i - ar{x})^k. ]
This uses the third central moment (m_3), so it is very sensitive to outliers.
def sample_skewness(x: np.ndarray) -> float:
"Moment-based sample skewness: m3 / m2^(3/2) using 1/n moments."
x = np.asarray(x, dtype=float).ravel()
x = x[np.isfinite(x)]
n = x.size
if n < 3:
return float("nan")
mean = x.mean()
centered = x - mean
m2 = np.mean(centered**2)
if m2 == 0:
return 0.0
m3 = np.mean(centered**3)
return float(m3 / (m2 ** 1.5))
n_demo = 4000
x_sym = rng.normal(0, 1, size=n_demo)
x_right = rng.exponential(scale=1.0, size=n_demo)
x_left = -rng.exponential(scale=1.0, size=n_demo)
skews = {
"Symmetric (normal)": sample_skewness(x_sym),
"Right-skewed (exponential)": sample_skewness(x_right),
"Left-skewed (-exponential)": sample_skewness(x_left),
}
skews
{'Symmetric (normal)': 0.003983194961544123,
'Right-skewed (exponential)': 1.979806008441165,
'Left-skewed (-exponential)': -1.901827954731168}
fig = make_subplots(
rows=1,
cols=3,
subplot_titles=[f"{k}<br>skew={v:.3f}" for k, v in skews.items()],
)
datasets = [
("Symmetric (normal)", x_sym),
("Right-skewed (exponential)", x_right),
("Left-skewed (-exponential)", x_left),
]
for col, (name, x) in enumerate(datasets, start=1):
fig.add_trace(
go.Histogram(
x=x,
nbinsx=80,
histnorm="probability density",
name=name,
marker=dict(line=dict(width=0)),
),
row=1,
col=col,
)
fig.update_layout(
height=320,
width=1000,
showlegend=False,
title_text="Skewness: symmetric vs right/left skew",
)
fig.update_xaxes(range=[-5, 5])
fig.show()
2) What the skewness test is (and what it is not)#
Goal#
We want to know whether the skewness we observe could plausibly be due to sampling noise when the underlying distribution is normal (which has population skewness 0).
Hypotheses#
For the standard two-sided skewness test:
(H_0): the population skewness is the same as the normal distribution (i.e. 0)
(H_1): the population skewness is not 0
You can also run one-sided versions:
alternative='greater': (H_1): skewness (> 0) (right-skew)alternative='less': (H_1): skewness (< 0) (left-skew)
What it is not#
It is not a complete normality test. A distribution can be symmetric (skew ≈ 0) but still very non-normal (e.g. heavy-tailed).
It does not tell you why data are skewed (mixtures, censoring, outliers, bounded support, etc.). You should always visualize the data.
3) D’Agostino’s skewness test: the core idea#
D’Agostino’s skewness test takes the sample skewness (b_2) and transforms it into a value (Z) that is approximately standard normal under (H_0):
[ Z pprox \mathcal{N}(0,1) \quad ext{(when the data are normal and } n\ge 8 ext{)}. ]
That lets us compute a p-value just like a z-test.
The transformation (as used by SciPy’s skewtest)#
Let (n) be the sample size and (b_2) the (moment-based) sample skewness. Define
[ y = b_2\sqrt{rac{(n+1)(n+3)}{6(n-2)}}. ]
Then compute
[ eta_2 = rac{3,(n^2 + 27n - 70)(n+1)(n+3)}{(n-2)(n+5)(n+7)(n+9)}, \quad W^2 = -1 + \sqrt{2(eta_2 - 1)}. ]
and
[ \delta = rac{1}{\sqrt{ frac{1}{2}\ln(W^2)}}, \quad lpha = \sqrt{rac{2}{W^2 - 1}}. ]
Finally,
[ Z = \delta,\operatorname{asinh}!\left(rac{y}{lpha} ight) = \delta,\ln\left(rac{y}{lpha} + \sqrt{\left(rac{y}{lpha} ight)^2 + 1} ight). ]
Why the sample size constraint?#
The approximation is derived under normality and is intended for (n\ge 8). For smaller samples, you should prefer resampling-based checks or simply treat skewness as an effect size.
4) NumPy-only implementation#
The functions below implement:
sample skewness (b_2)
D’Agostino’s transformed statistic (Z)
p-values using a NumPy approximation to the normal CDF (via an
erfapproximation)
The goal here is transparency, not squeezing out every last bit of numerical accuracy.
def erf_approx(x: np.ndarray) -> np.ndarray:
"Vectorized erf approximation (Abramowitz & Stegun 7.1.26)."
x = np.asarray(x, dtype=float)
sign = np.sign(x)
ax = np.abs(x)
p = 0.3275911
t = 1.0 / (1.0 + p * ax)
a1 = 0.254829592
a2 = -0.284496736
a3 = 1.421413741
a4 = -1.453152027
a5 = 1.061405429
poly = (((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t)
y = 1.0 - poly * np.exp(-ax * ax)
return sign * y
def norm_cdf(z: np.ndarray) -> np.ndarray:
"Standard normal CDF Φ(z) using erf approximation."
z = np.asarray(z, dtype=float)
return 0.5 * (1.0 + erf_approx(z / np.sqrt(2.0)))
def pvalue_from_z(z: np.ndarray, alternative: str = "two-sided") -> np.ndarray:
"Convert a z-statistic to a p-value under N(0,1)."
z = np.asarray(z, dtype=float)
if alternative == "two-sided":
p = 2.0 * norm_cdf(-np.abs(z))
elif alternative == "greater":
p = norm_cdf(-z)
elif alternative == "less":
p = norm_cdf(z)
else:
raise ValueError("alternative must be one of: 'two-sided', 'greater', 'less'")
return np.clip(p, 0.0, 1.0)
def dagostino_z_from_skewness(b2: np.ndarray, n: int) -> np.ndarray:
"D'Agostino transform from sample skewness b2 to Z (approx N(0,1) under H0)."
if n < 8:
raise ValueError("D'Agostino skewness test requires n >= 8")
b2 = np.asarray(b2, dtype=float)
nf = float(n)
y = b2 * np.sqrt(((nf + 1.0) * (nf + 3.0)) / (6.0 * (nf - 2.0)))
beta2 = (
3.0
* (nf**2 + 27.0 * nf - 70.0)
* (nf + 1.0)
* (nf + 3.0)
/ ((nf - 2.0) * (nf + 5.0) * (nf + 7.0) * (nf + 9.0))
)
w2 = -1.0 + np.sqrt(2.0 * (beta2 - 1.0))
delta = 1.0 / np.sqrt(0.5 * np.log(w2))
alpha = np.sqrt(2.0 / (w2 - 1.0))
# asinh(u) == log(u + sqrt(u^2 + 1)) but is more numerically stable
return delta * np.arcsinh(y / alpha)
def dagostino_skewtest(x: np.ndarray, alternative: str = "two-sided") -> dict:
"D'Agostino skewness test (NumPy-only implementation)."
x = np.asarray(x, dtype=float).ravel()
x = x[np.isfinite(x)]
n = x.size
if n < 8:
raise ValueError("D'Agostino skewness test requires at least 8 observations")
b2 = sample_skewness(x)
z = dagostino_z_from_skewness(b2, n)
p = pvalue_from_z(z, alternative=alternative)
return {
"n": int(n),
"skewness": float(b2),
"z": float(np.asarray(z)),
"pvalue": float(np.asarray(p)),
"alternative": alternative,
}
5) Worked examples + interpretation#
Interpretation guide (two-sided):
Small p-value (e.g. < 0.05) → evidence the data are asymmetric (skewness ≠ 0) relative to normal sampling noise.
The sign of Z tells you the direction:
Z > 0→ right-skew (long right tail)Z < 0→ left-skew (long left tail)
Remember: with large (n), even tiny skewness can become “statistically significant”. Always look at the effect size (the skewness value) and the plot.
def summarize_test(x: np.ndarray, name: str, alternative: str = "two-sided") -> dict:
res = dagostino_skewtest(x, alternative=alternative)
return {
"dataset": name,
"n": res["n"],
"skewness": res["skewness"],
"z": res["z"],
"pvalue": res["pvalue"],
"alternative": res["alternative"],
}
n = 60
x1 = rng.normal(0, 1, size=n)
x2 = rng.exponential(1.0, size=n)
x3 = -rng.exponential(1.0, size=n)
results = [
summarize_test(x1, "Normal sample"),
summarize_test(x2, "Right-skew (exponential)"),
summarize_test(x3, "Left-skew (-exponential)"),
]
results
[{'dataset': 'Normal sample',
'n': 60,
'skewness': 0.33403897478489747,
'z': 1.1391470464804974,
'pvalue': 0.25464193619392184,
'alternative': 'two-sided'},
{'dataset': 'Right-skew (exponential)',
'n': 60,
'skewness': 1.0681651963774472,
'z': 3.204306467065742,
'pvalue': 0.001354010697767749,
'alternative': 'two-sided'},
{'dataset': 'Left-skew (-exponential)',
'n': 60,
'skewness': -2.314018629424382,
'z': -5.399094653879995,
'pvalue': 6.713003297686981e-08,
'alternative': 'two-sided'}]
fig = make_subplots(
rows=1,
cols=3,
subplot_titles=[
f"{r['dataset']}<br>skew={r['skewness']:.3f}, Z={r['z']:.2f}, p={r['pvalue']:.3g}" for r in results
],
)
for col, (x, r) in enumerate([(x1, results[0]), (x2, results[1]), (x3, results[2])], start=1):
fig.add_trace(
go.Histogram(
x=x,
nbinsx=30,
histnorm="probability density",
name=r["dataset"],
marker=dict(line=dict(width=0)),
),
row=1,
col=col,
)
fig.update_layout(
height=320,
width=1100,
showlegend=False,
title_text="Skewness test on small samples (n=60)",
)
fig.show()
One-sided alternatives#
If you have a directional question, you can use a one-sided test:
greater: “is it right-skewed?” (skewness > 0)less: “is it left-skewed?” (skewness < 0)
One-sided p-values are smaller only when the observed skew is in the hypothesized direction.
x = rng.exponential(1.0, size=60)
res_two = dagostino_skewtest(x, alternative="two-sided")
res_greater = dagostino_skewtest(x, alternative="greater")
res_less = dagostino_skewtest(x, alternative="less")
{
"two-sided": res_two,
"greater": res_greater,
"less": res_less,
}
{'two-sided': {'n': 60,
'skewness': 2.616686399904669,
'z': 5.78336467838644,
'pvalue': 7.344590824409636e-09,
'alternative': 'two-sided'},
'greater': {'n': 60,
'skewness': 2.616686399904669,
'z': 5.78336467838644,
'pvalue': 3.672295412204818e-09,
'alternative': 'greater'},
'less': {'n': 60,
'skewness': 2.616686399904669,
'z': 5.78336467838644,
'pvalue': 0.9999999963277046,
'alternative': 'less'}}
6) Why sample size matters (p-value vs skewness)#
Because the test statistic is a standardized measure of skewness, larger samples make it easier to detect small asymmetries.
The plot below shows the two-sided p-value as a function of sample skewness (b_2) for different (n).
b2_grid = np.linspace(-2.5, 2.5, 501)
fig = go.Figure()
for n in [10, 30, 100, 500]:
z_grid = dagostino_z_from_skewness(b2_grid, n=n)
p_grid = pvalue_from_z(z_grid, alternative="two-sided")
fig.add_trace(go.Scatter(x=b2_grid, y=p_grid, mode="lines", name=f"n={n}"))
fig.add_hline(y=0.05, line_dash="dash", line_color="black")
fig.update_layout(
title="Two-sided p-value vs sample skewness (D'Agostino transform)",
xaxis_title="sample skewness b2",
yaxis_title="p-value",
yaxis_type="log",
height=420,
width=900,
)
fig.show()
7) Does Z really look standard normal under H0?#
Under the null (normal data), the transformed statistic (Z) is designed to be approximately (\mathcal{N}(0,1)).
We can sanity-check that with a small Monte Carlo simulation.
def skewness_vectorized(samples: np.ndarray) -> np.ndarray:
"Vectorized skewness for samples shaped (n_sims, n)."
samples = np.asarray(samples, dtype=float)
mean = samples.mean(axis=1, keepdims=True)
centered = samples - mean
m2 = np.mean(centered**2, axis=1)
m3 = np.mean(centered**3, axis=1)
out = np.zeros_like(m2)
mask = m2 > 0
out[mask] = m3[mask] / (m2[mask] ** 1.5)
out[~mask] = 0.0
return out
n = 20
n_sims = 10_000
samples = rng.normal(0, 1, size=(n_sims, n))
b2 = skewness_vectorized(samples)
z = dagostino_z_from_skewness(b2, n=n)
z_mean = float(np.mean(z))
z_std = float(np.std(z, ddof=0))
rejection_rate = float(np.mean(pvalue_from_z(z, alternative="two-sided") < 0.05))
{"z_mean": z_mean, "z_std": z_std, "rejection_rate@0.05": rejection_rate}
{'z_mean': 0.010093399552716806,
'z_std': 0.9999562884320329,
'rejection_rate@0.05': 0.0476}
# Histogram of Z with a standard normal PDF overlay
x_grid = np.linspace(-4, 4, 400)
pdf = (1.0 / np.sqrt(2 * np.pi)) * np.exp(-0.5 * x_grid**2)
fig = go.Figure()
fig.add_trace(
go.Histogram(
x=z,
nbinsx=60,
histnorm="probability density",
name="Simulated Z (H0: normal)",
opacity=0.75,
marker=dict(line=dict(width=0)),
)
)
fig.add_trace(go.Scatter(x=x_grid, y=pdf, mode="lines", name="N(0,1) PDF", line=dict(color="black")))
fig.update_layout(
title=f"Under H0, Z is approximately N(0,1) (n={n}, sims={n_sims})",
xaxis_title="Z",
yaxis_title="density",
height=420,
width=900,
)
fig.show()
8) Practical notes and pitfalls#
Outliers can dominate: skewness uses ((x-ar{x})^3). One extreme point can flip the conclusion.
Statistical vs practical significance: with large (n), tiny skewness can yield tiny p-values.
Not a full normality test: skewness ≈ 0 does not imply normality.
Independence matters: if observations are dependent (time series, clustered data), p-values can be misleading.
Use visuals: combine the test with histograms/QQ-plots and domain knowledge.
9) Optional: compare with SciPy#
If you have SciPy available, you can verify that the NumPy implementation matches scipy.stats.skewtest.
x = rng.lognormal(mean=0.0, sigma=0.8, size=200)
ours = dagostino_skewtest(x)
try:
from scipy.stats import skewtest
scipy_res = skewtest(x)
comparison = {
"ours": ours,
"scipy": {"z": float(scipy_res.statistic), "pvalue": float(scipy_res.pvalue)},
"abs_diff_z": float(abs(ours["z"] - scipy_res.statistic)),
"abs_diff_p": float(abs(ours["pvalue"] - scipy_res.pvalue)),
}
except Exception as e:
comparison = {"ours": ours, "scipy": None, "error": repr(e)}
comparison
{'ours': {'n': 200,
'skewness': 2.5556762139133857,
'z': 9.304601238442498,
'pvalue': 0.0,
'alternative': 'two-sided'},
'scipy': {'z': 9.304601238442498, 'pvalue': 1.3449606122843268e-20},
'abs_diff_z': 0.0,
'abs_diff_p': 1.3449606122843268e-20}
Exercises#
Simulate data from a symmetric but heavy-tailed distribution (e.g. Student t) and see how often the skewness test rejects.
Create a mixture of two normals with different means and observe how skewness changes.
For a fixed skewness, increase (n) and see how quickly the p-value shrinks.
References#
R. B. D’Agostino et al. (1990), A suggestion for using powerful and informative tests of normality
SciPy documentation for
scipy.stats.skewtest